Stepwise Multiple Regressions
Leah Thornton Hampton
2/1/2021
1 Set Up
#Custom function that finds factors with only 1 level
#Source: https://stackoverflow.com/questions/44200195/how-to-debug-contrasts-can-be-applied-only-to-factors-with-2-or-more-levels-er
debug_contr_error <- function (dat, subset_vec = NULL) {
if (!is.null(subset_vec)) {
## step 0
if (mode(subset_vec) == "logical") {
if (length(subset_vec) != nrow(dat)) {
stop("'logical' `subset_vec` provided but length does not match `nrow(dat)`")
}
subset_log_vec <- subset_vec
} else if (mode(subset_vec) == "numeric") {
## check range
ran <- range(subset_vec)
if (ran[1] < 1 || ran[2] > nrow(dat)) {
stop("'numeric' `subset_vec` provided but values are out of bound")
} else {
subset_log_vec <- logical(nrow(dat))
subset_log_vec[as.integer(subset_vec)] <- TRUE
}
} else {
stop("`subset_vec` must be either 'logical' or 'numeric'")
}
dat <- base::subset(dat, subset = subset_log_vec)
} else {
## step 1
dat <- stats::na.omit(dat)
}
if (nrow(dat) == 0L) warning("no complete cases")
## step 2
var_mode <- sapply(dat, mode)
if (any(var_mode %in% c("complex", "raw"))) stop("complex or raw not allowed!")
var_class <- sapply(dat, class)
if (any(var_mode[var_class == "AsIs"] %in% c("logical", "character"))) {
stop("matrix variables with 'AsIs' class must be 'numeric'")
}
ind1 <- which(var_mode %in% c("logical", "character"))
dat[ind1] <- lapply(dat[ind1], as.factor)
## step 3
fctr <- which(sapply(dat, is.factor))
if (length(fctr) == 0L) warning("no factor variables to summary")
ind2 <- if (length(ind1) > 0L) fctr[-ind1] else fctr
dat[ind2] <- lapply(dat[ind2], base::droplevels.factor)
## step 4
lev <- lapply(dat[fctr], base::levels.default)
nl <- lengths(lev)
## return
list(nlevels = nl, levels = lev)
}#Variables to put into model
- Dose (counts)
- Dose (mass)
- Dose (volume)
- Size (continuous)
- Size (binned)
- Polymer type
- Shape
- Organism Group (this is being used as a temporary substitute for body size)
- Life Stage
- Level of Biological Organization
- Exposure Duration
- Acute/Chronic (this values are only complete for fish, molluscs, crustaceans and algae)
- Charge (negative/positive - categorical)
- Zeta potential
- Particle Source (commercial, generated in-house, or mmodified commercial particles e.g., milling)
Response Variables * Effect (Y/N) - Binary Data * Effect (Y/N) - Binary Data x Effect Score (binned organism level effects)
#Continous Variable Distributions
There are several continuous variables that we want to feed into our model. Before doing so, we’re going to check the distribution to see if any of them are skewed and need to be transformed.
Due to skewed data, the following categories are log10 transformed before modeling:
- Dose (counts)
- Dose (mass)
- Dose (volume)
- Size (continuous)
- Exposure Duration
2 Model: Kitchen Sink
All Independent Variables, Response Variable: Effect (Y/N)
Full Model
## $nlevels
## shape_f org_f life_f bio_f acute.chronic_f
## 3 3 3 5 2
##
## $levels
## $levels$shape_f
## [1] "Fiber" "Fragment" "Sphere"
##
## $levels$org_f
## [1] "Crustacea" "Fish" "Mollusca"
##
## $levels$life_f
## [1] "Adult" "Early" "Juvenile"
##
## $levels$bio_f
## [1] "Cell" "Organism" "Population" "Subcell" "Tissue"
##
## $levels$acute.chronic_f
## [1] "Acute" "Chronic"##
## Call:
## lm(formula = effect_10 ~ ., data = aoc_setup_select_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6552 -0.3369 -0.1877 0.4716 1.1018
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.051687 0.150037 7.010 2.87e-12
## log.dose.particles.mL.master -0.112899 0.125461 -0.900 0.36825
## log.dose.mg.L.master 0.051598 0.008582 6.012 2.03e-09
## log.dose.um3.mL.master 0.093169 0.125761 0.741 0.45884
## log.size.length.um.used.for.conversions -0.361561 0.376854 -0.959 0.33742
## shape_fFragment -0.261101 0.046653 -5.597 2.36e-08
## shape_fSphere -0.271240 0.101701 -2.667 0.00769
## org_fFish -0.071743 0.023682 -3.029 0.00247
## org_fMollusca -0.165614 0.026550 -6.238 4.98e-10
## life_fEarly -0.054630 0.022028 -2.480 0.01319
## life_fJuvenile -0.110397 0.023656 -4.667 3.18e-06
## bio_fOrganism -0.341509 0.047869 -7.134 1.18e-12
## bio_fPopulation -0.616304 0.095914 -6.426 1.50e-10
## bio_fSubcell -0.124617 0.045353 -2.748 0.00603
## bio_fTissue -0.167895 0.055423 -3.029 0.00247
## log.exposure.duration.d 0.058345 0.014896 3.917 9.15e-05
## acute.chronic_fChronic 0.051101 0.023834 2.144 0.03210
##
## (Intercept) ***
## log.dose.particles.mL.master
## log.dose.mg.L.master ***
## log.dose.um3.mL.master
## log.size.length.um.used.for.conversions
## shape_fFragment ***
## shape_fSphere **
## org_fFish **
## org_fMollusca ***
## life_fEarly *
## life_fJuvenile ***
## bio_fOrganism ***
## bio_fPopulation ***
## bio_fSubcell **
## bio_fTissue **
## log.exposure.duration.d ***
## acute.chronic_fChronic *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4382 on 3385 degrees of freedom
## Multiple R-squared: 0.1104, Adjusted R-squared: 0.1062
## F-statistic: 26.25 on 16 and 3385 DF, p-value: < 2.2e-16Stepwise Model - Both Directions
##
## Call:
## lm(formula = effect_10 ~ log.dose.particles.mL.master + log.dose.mg.L.master +
## log.size.length.um.used.for.conversions + shape_f + org_f +
## life_f + bio_f + log.exposure.duration.d + acute.chronic_f,
## data = aoc_setup_select_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6616 -0.3384 -0.1883 0.4722 1.1004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.961328 0.087373 11.003 < 2e-16
## log.dose.particles.mL.master -0.020160 0.008385 -2.404 0.01626
## log.dose.mg.L.master 0.051909 0.008571 6.056 1.55e-09
## log.size.length.um.used.for.conversions -0.083029 0.025820 -3.216 0.00131
## shape_fFragment -0.268834 0.045467 -5.913 3.70e-09
## shape_fSphere -0.205081 0.048657 -4.215 2.57e-05
## org_fFish -0.071653 0.023680 -3.026 0.00250
## org_fMollusca -0.164836 0.026527 -6.214 5.80e-10
## life_fEarly -0.053066 0.021925 -2.420 0.01556
## life_fJuvenile -0.109698 0.023636 -4.641 3.59e-06
## bio_fOrganism -0.341726 0.047865 -7.139 1.14e-12
## bio_fPopulation -0.616644 0.095906 -6.430 1.46e-10
## bio_fSubcell -0.124734 0.045350 -2.750 0.00598
## bio_fTissue -0.168007 0.055420 -3.032 0.00245
## log.exposure.duration.d 0.058166 0.014893 3.905 9.59e-05
## acute.chronic_fChronic 0.050621 0.023824 2.125 0.03368
##
## (Intercept) ***
## log.dose.particles.mL.master *
## log.dose.mg.L.master ***
## log.size.length.um.used.for.conversions **
## shape_fFragment ***
## shape_fSphere ***
## org_fFish **
## org_fMollusca ***
## life_fEarly *
## life_fJuvenile ***
## bio_fOrganism ***
## bio_fPopulation ***
## bio_fSubcell **
## bio_fTissue **
## log.exposure.duration.d ***
## acute.chronic_fChronic *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4381 on 3386 degrees of freedom
## Multiple R-squared: 0.1102, Adjusted R-squared: 0.1063
## F-statistic: 27.96 on 15 and 3386 DF, p-value: < 2.2e-163 Model: Crustacean Fitness, Organism Level Endpoints Only
3.1 Dose (mass)
Full Model
##
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8447 -0.7063 -0.6178 -0.5142 2.0485
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -1.425e+00 7.388e-02
## logdose.mg.L.master 1.807e-01 4.837e-02
## size.length.um.used.for.conversions -3.807e-04 1.785e-04
## logdose.mg.L.master:size.length.um.used.for.conversions 1.057e-04 4.555e-05
## z value Pr(>|z|)
## (Intercept) -19.288 < 2e-16 ***
## logdose.mg.L.master 3.737 0.000187 ***
## size.length.um.used.for.conversions -2.133 0.032944 *
## logdose.mg.L.master:size.length.um.used.for.conversions 2.320 0.020367 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1272.9 on 1307 degrees of freedom
## AIC: 1280.9
##
## Number of Fisher Scoring iterations: 4#Plots
ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)3.2 Dose (volume)
Full Model
##
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8116 -0.6866 -0.6282 -0.5292 2.0869
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -2.176e+00 2.574e-01
## logdose.um3.mL.master 1.262e-01 4.049e-02
## size.length.um.used.for.conversions -1.176e-03 4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions 1.259e-04 4.533e-05
## z value Pr(>|z|)
## (Intercept) -8.451 < 2e-16 ***
## logdose.um3.mL.master 3.118 0.00182 **
## size.length.um.used.for.conversions -2.635 0.00842 **
## logdose.um3.mL.master:size.length.um.used.for.conversions 2.777 0.00549 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1277.6 on 1307 degrees of freedom
## AIC: 1285.6
##
## Number of Fisher Scoring iterations: 4#Plots
ggPredict(m1_crust_model_dose,
colorn=200,
# jitter=TRUE,
point = FALSE,
interactive = TRUE,
show.summary = TRUE,
#digits = 3,
se = FALSE)##
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8116 -0.6866 -0.6282 -0.5292 2.0869
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -2.176e+00 2.574e-01
## logdose.um3.mL.master 1.262e-01 4.049e-02
## size.length.um.used.for.conversions -1.176e-03 4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions 1.259e-04 4.533e-05
## z value Pr(>|z|)
## (Intercept) -8.451 < 2e-16 ***
## logdose.um3.mL.master 3.118 0.00182 **
## size.length.um.used.for.conversions -2.635 0.00842 **
## logdose.um3.mL.master:size.length.um.used.for.conversions 2.777 0.00549 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1277.6 on 1307 degrees of freedom
## AIC: 1285.6
##
## Number of Fisher Scoring iterations: 4#particle couunt
#Select the columns that we want to feed into the model for simplicity
m1_crust_part <- aoc_setup %>%
filter(!effect_f == "NA") %>%
#filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>%
#filter(organism.group == "Crustacea") %>%
filter(!acute.chronic_f == "Unavailable") %>%
filter(bio.org == "organism") %>%
filter(!environment == "Terrestrial") %>%
filter(lvl1_f == "Fitness") %>%
filter(!size_f == "Not Reported") %>%
mutate(logdose.particles.mL.master = log10(dose.particles.mL.master)) %>%
# filter(acute.chronic_f == "Acute") %>%
dplyr::select(c(effect_10, logdose.particles.mL.master, size.length.um.used.for.conversions)) %>%
drop_na()
#Fit the full model
m1_crust_model_particles <- glm(effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions, data = m1_crust_part, na.action = "na.exclude", family = "binomial")
summary(m1_crust_model_particles)##
## Call:
## glm(formula = effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust_part, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5942 -0.6724 -0.6214 -0.5551 2.0530
##
## Coefficients:
## Estimate
## (Intercept) -1.808e+00
## logdose.particles.mL.master 7.178e-02
## size.length.um.used.for.conversions 2.881e-04
## logdose.particles.mL.master:size.length.um.used.for.conversions 1.331e-04
## Std. Error
## (Intercept) 1.355e-01
## logdose.particles.mL.master 2.177e-02
## size.length.um.used.for.conversions 6.364e-05
## logdose.particles.mL.master:size.length.um.used.for.conversions 4.499e-05
## z value
## (Intercept) -13.342
## logdose.particles.mL.master 3.297
## size.length.um.used.for.conversions 4.528
## logdose.particles.mL.master:size.length.um.used.for.conversions 2.959
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## logdose.particles.mL.master 0.000978 ***
## size.length.um.used.for.conversions 5.96e-06 ***
## logdose.particles.mL.master:size.length.um.used.for.conversions 0.003090 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1280.0 on 1307 degrees of freedom
## AIC: 1288
##
## Number of Fisher Scoring iterations: 4plot
ggPredict(m1_crust_model_particles, colorn=1000, interactive = TRUE)3.2.1 volume
#particle couunt
#Select the columns that we want to feed into the model for simplicity
m1_crust_volume <- aoc_setup %>%
filter(!effect_f == "NA") %>%
#filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>%
#filter(organism.group == "Crustacea") %>%
filter(!acute.chronic_f == "Unavailable") %>%
filter(bio.org == "organism") %>%
filter(!environment == "Terrestrial") %>%
filter(lvl1_f == "Fitness") %>%
filter(!size_f == "Not Reported") %>%
mutate(logdose.um3.mL.master = log10(dose.um3.mL.master)) %>%
# filter(acute.chronic_f == "Acute") %>%
dplyr::select(c(effect_10, logdose.um3.mL.master, size.length.um.used.for.conversions)) %>%
drop_na()
#Fit the full model
m1_crust_model_volume <- glm(effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions, data = m1_crust_volume, na.action = "na.exclude", family = "binomial")
summary(m1_crust_model_volume)##
## Call:
## glm(formula = effect_10 ~ logdose.um3.mL.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust_volume, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8116 -0.6866 -0.6282 -0.5292 2.0869
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -2.176e+00 2.574e-01
## logdose.um3.mL.master 1.262e-01 4.049e-02
## size.length.um.used.for.conversions -1.176e-03 4.464e-04
## logdose.um3.mL.master:size.length.um.used.for.conversions 1.259e-04 4.533e-05
## z value Pr(>|z|)
## (Intercept) -8.451 < 2e-16 ***
## logdose.um3.mL.master 3.118 0.00182 **
## size.length.um.used.for.conversions -2.635 0.00842 **
## logdose.um3.mL.master:size.length.um.used.for.conversions 2.777 0.00549 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1277.6 on 1307 degrees of freedom
## AIC: 1285.6
##
## Number of Fisher Scoring iterations: 4plot
ggPredict(m1_crust_model_volume, colorn=500, interactive = TRUE)#alternative method using ggplot
#get probabilities
m1 <- m1_crust %>%
mutate(prob = predict(m1_crust_model_dose,
type = "response"))
# plot
ggplot(m1, aes(x = logdose.mg.L.master, y = effect_10, color = size.length.um.used.for.conversions))+
geom_jitter() +
geom_smooth(method = 'glm',
aes(x = logdose.mg.L.master, y = prob), #use predicted NOEC probability
method.args = list(family = binomial(link = 'logit')),
se = FALSE, color = 'red') #+
# geom_smooth(method = 'glm',
# method.args = list(family = binomial(link = 'probit')),
# se = FALSE, color = 'green') +
# stat_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE)ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)3.3 Dose (count)
Full Model
##
## Call:
## glm(formula = effect_10 ~ logdose.particles.mL.master * size.length.um.used.for.conversions,
## family = "binomial", data = m1_crust, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5942 -0.6724 -0.6214 -0.5551 2.0530
##
## Coefficients:
## Estimate
## (Intercept) -1.808e+00
## logdose.particles.mL.master 7.178e-02
## size.length.um.used.for.conversions 2.881e-04
## logdose.particles.mL.master:size.length.um.used.for.conversions 1.331e-04
## Std. Error
## (Intercept) 1.355e-01
## logdose.particles.mL.master 2.177e-02
## size.length.um.used.for.conversions 6.364e-05
## logdose.particles.mL.master:size.length.um.used.for.conversions 4.499e-05
## z value
## (Intercept) -13.342
## logdose.particles.mL.master 3.297
## size.length.um.used.for.conversions 4.528
## logdose.particles.mL.master:size.length.um.used.for.conversions 2.959
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## logdose.particles.mL.master 0.000978 ***
## size.length.um.used.for.conversions 5.96e-06 ***
## logdose.particles.mL.master:size.length.um.used.for.conversions 0.003090 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1311.5 on 1310 degrees of freedom
## Residual deviance: 1280.0 on 1307 degrees of freedom
## AIC: 1288
##
## Number of Fisher Scoring iterations: 4require(ggiraphExtra)
#Plots
ggPredict(m1_crust_model_dose,colorn=100,jitter=FALSE, interactive = TRUE)3.4 Discrete predictor and dose
##
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size_f, family = "binomial",
## data = m1_crust, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7085 -0.7062 -0.6167 -0.4255 2.3881
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.57691 0.21094 -7.476 7.67e-14 ***
## logdose.mg.L.master 0.40501 0.16344 2.478 0.0132 *
## size_f100nm < 1µm -0.01864 0.37454 -0.050 0.9603
## size_f1µm < 100µm 0.17395 0.22858 0.761 0.4467
## size_f100µm < 1mm -0.25876 0.35626 -0.726 0.4676
## size_f1mm < 5mm -0.87511 0.64223 -1.363 0.1730
## logdose.mg.L.master:size_f100nm < 1µm 0.45050 0.31359 1.437 0.1508
## logdose.mg.L.master:size_f1µm < 100µm -0.24482 0.17319 -1.414 0.1575
## logdose.mg.L.master:size_f100µm < 1mm -0.23277 0.22994 -1.012 0.3114
## logdose.mg.L.master:size_f1mm < 5mm 0.04129 0.21904 0.189 0.8505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1358.5 on 1347 degrees of freedom
## Residual deviance: 1303.8 on 1338 degrees of freedom
## AIC: 1323.8
##
## Number of Fisher Scoring iterations: 5Plotted below.
discrete <- ggpredict(m1_crust_model_discrete,terms = c("logdose.mg.L.master", "size_f"))
plot(discrete, facet = TRUE, add.data = TRUE, colors = "quadro", alpha = 0.2, dodge = 0.2)
#title = "All Aquatic Organisms")3.4.0.1 Acute only crustacea
##
## Call:
## glm(formula = effect_10 ~ logdose.mg.L.master * size_f + logdose.particles.mL.master *
## size_f, family = "binomial", data = m1_crust_acute, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4072 -0.8096 -0.6639 -0.2631 2.1764
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -30.504 9.991 -3.053
## logdose.mg.L.master -2.914 1.075 -2.711
## size_f100nm < 1µm 29.153 10.270 2.839
## size_f1µm < 100µm 28.853 9.996 2.886
## size_f100µm < 1mm 28.939 10.003 2.893
## logdose.particles.mL.master 2.912 1.001 2.908
## logdose.mg.L.master:size_f100nm < 1µm 3.637 1.165 3.123
## logdose.mg.L.master:size_f1µm < 100µm 3.002 1.078 2.785
## logdose.mg.L.master:size_f100µm < 1mm 2.298 1.148 2.002
## size_f100nm < 1µm:logdose.particles.mL.master -2.905 1.047 -2.773
## size_f1µm < 100µm:logdose.particles.mL.master -2.767 1.004 -2.756
## size_f100µm < 1mm:logdose.particles.mL.master -1.988 1.153 -1.725
## Pr(>|z|)
## (Intercept) 0.00227 **
## logdose.mg.L.master 0.00672 **
## size_f100nm < 1µm 0.00453 **
## size_f1µm < 100µm 0.00390 **
## size_f100µm < 1mm 0.00382 **
## logdose.particles.mL.master 0.00364 **
## logdose.mg.L.master:size_f100nm < 1µm 0.00179 **
## logdose.mg.L.master:size_f1µm < 100µm 0.00535 **
## logdose.mg.L.master:size_f100µm < 1mm 0.04526 *
## size_f100nm < 1µm:logdose.particles.mL.master 0.00555 **
## size_f1µm < 100µm:logdose.particles.mL.master 0.00585 **
## size_f100µm < 1mm:logdose.particles.mL.master 0.08459 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 945.92 on 853 degrees of freedom
## Residual deviance: 904.63 on 842 degrees of freedom
## AIC: 928.63
##
## Number of Fisher Scoring iterations: 5The above glm is visualized below.
acute_discrete <- ggpredict(m1_crust_acute_model_discrete,terms = c(#"logdose.mg.L.master",
"logdose.particles.mL.master",
"size_f")) #ensure discrete factor is second
plot(acute_discrete)#, facet = TRUE, add.data = TRUE)
3.4.0.2 Alt Method using GG Plot
This approach achieves a similiar product as using the ggPredict() function, except it relies on ggplot(), which is more malleable and transparent. The general steps are to first create a new dataframe over 1000 values of size using expand.grid() then use predict() and plot() with geom_line() and colour=size.
require(viridis)
#filtered dataset
#m1_crust
#model
#summary(m1_crust_model_dose )
#generate distribution of data
mockData <- expand.grid(size.length.um.used.for.conversions = seq(0.034, 5000, 1),
logdose.particles.mL.master = seq(-4.195, 12.650,0.1))
mockData$effect_10 <- predict(m1_crust_model_dose,
mockData,
type = "response")
#plot
mockData %>%
ggplot(aes(x = logdose.particles.mL.master,
y = effect_10,
color = size.length.um.used.for.conversions)) +
geom_line() +
scale_color_gradientn(colors = topo.colors(7)) +
#scale_color_viridis_c(option = "A") +
dark_theme_bw()
3.5 Survival package
3.5.1 Cumulative Hazard by Time
#include time
m1_crust <- aoc_setup %>%
filter(!effect_f == "NA") %>%
#filter(effect.metric == c("HONEC", "LOEC", "NOEC")) %>%
#filter(organism.group == "Crustacea") %>%
filter(!acute.chronic_f == "Unavailable") %>%
filter(bio.org == "organism") %>%
filter(!environment == "Terrestrial") %>%
filter(lvl1_f == "Fitness") %>%
filter(!size_f == "Not Reported") %>%
mutate(logdose.mg.L.master = log10(dose.mg.L.master)) %>%
# filter(acute.chronic_f == "Acute") %>%
dplyr::select(c(effect_10, logdose.mg.L.master, size_f, exposure.duration.d)) %>%
drop_na()
require(survival)
survival <- survival::coxph(survival::Surv(exposure.duration.d, effect_10) ~ logdose.mg.L.master + size_f, data = m1_crust)
#cumulative hazard
cumHaz <- ggpredict(survival,terms = c("logdose.mg.L.master", "size_f"), type = "cumhaz")
plot(cumHaz, facet = TRUE, add.data = TRUE, colors = "flat")
### Probability of Survival by Time
#survival probability over time
pr <- ggpredict(survival, c("logdose.mg.L.master", "size_f"), type = "surv")
plot(pr, colors = "social")